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Abstract This paper develops a polynomial normal transformation model, whereby various non¬ 
normal probability distributions can be simulated by the standard normal distribution. Two meth¬ 
ods are presented to determine the coefficients of polynomial model: (1) probability weighted mo¬ 
ment (PWM) matching (2) percentile matching. Compared to the existing raw moment or L-moment 
matching, the proposed methods are more computationally convenient, and can be used to estimate 
the coefficients of polynomial model with a higher degree. Furthermore, for two correlated random 
variables, a polynomial equation is derived to estimate the equivalent correlation coefficient in 
standard normal space, and random vector with non-normal marginal distributions and prescribed 
correlation matrix can be generated. Finally, numerical examples are worked to demonstrate the 
proposed method. 

Keywords Polynomial normal transformation • NORTA algorithm ■ Probability weighted 
moment • Percentile • Correlation coefficient 


1 Introduction 

In order to generate non-normal random variables, many families of distributions have been pro¬ 
posed. [Tadikamalla(1980)] summarizes the related works, and compares six procedures in terms of 
the efficiency, simplicity and generality. Among these approaches, two standard normal distribution 
based transformation models: Johnson system and third-order polynomial normal transformation 
(TPNT) technique, demonstrate the potential utility for generating correlated random variables. 

[Johnson(1949)] discusses the feasibility of simulating continuous non-normal distributions by 
standard normal distribution, and develops a four-parameter transformation system. The John¬ 
son system comprises four transformation models: Sn (Normal distribution), Sl (Lognormal fam¬ 
ily), Su (Unbounded family) and Sb (Bounded family), which can be employed to approximate 
many non-normal distributions. For the detailed introduction to Johnson system, one can refer to 
literature[DeBrota et al(1989)]. 


Qing Xiao 

School of Mechatronic Engineering and Automation, Shanghai University, Shanghai, 200072, China 
E-mail: xaoshaoying@shu.edu.cn 





2 


Qing Xiao 


Johnson system has been widely used to generate non-normal random numbers, but it shows 
inadequacy of accommodating the dependence among random variables. If samples of the correlated 
random vector are available, the Johnson system is workable. By normalizing each sample value 
through the corresponding transformation model, the correlation matrix in the normal space can 
be calculated directly ( see Section 3.2 of reference [Stanfield et al(1996)]. However, this approach 
cannot be used to generate random vector with assigned correlation matrix. To achieve this goal, 
it is necessary to analytically estimate the correlation coefficient p z in normal space for a given 
correlation coefficient p x between two random variables. 


Compared with Johnson system, the dependence among random variables can be easily handled 
by TPNT technique. For a given p x , the analytical equation can be derived to estimate p 2 , and ran¬ 
dom vector with desired correlation matrix can be conveniently generated [Vale and Maurelli(1983), 
Headick and Sanwilowsky(1999)]. The TPNT is firstly introduced by [Fleishman(1978)], this tech¬ 
nique generates non-normal deviates by a polynomial transformation of standard normal ones. A 
detailed investigation of TPNT is given in reference[Chen and Tung(2003)], it is found there are 
many probability distributions that TPNT technique cannot simulate. 


To expand the generality of the polynomial transformation technique, a fifth-order polynomial 
normal transformation (FPNT) technique is developed in reference [Headick(2002)]. Due to the 
ability to control the extra fifth and sixth order moments, FPNT can approximate some distributions 
that are difficult for TPNT[Headrick and Kowalchuk(2007)], but it requires tedious mathematical 
computation to determine the coefficients of FPNT. By equating the first six raw moments of FPNT 
model with those of the modeled distribution, a system of nonlinear equations is established. The 
coefficients are estimated by solving the equations numerically. In pursuit of a simpler parameter 
estimation method, L-moments are employed to characterize FPNT model[Headick(2011)], and the 
analytical expression of coefficients are obtained. However, because of the tedious mathematical 
derivation, this L-moments based method is unable to estimate the coefficients of a polynomial 
model with a higher degree, which may be used to simulate more non-normal distributions. 


In this paper, the polynomial normal transformation model is extended to a much higher degree. 
Two methods: probability weighted moment (PWM) matching and percentile matching, are devel¬ 
oped to assess the coefficients of polynomial model. By equating PWMs of the polynomial model 
with those of the modeled distribution, a system of linear equations solved for the coefficients is 
established. For the percentile matching, the polynomial model can be established by a least square 
method, and a numerical investigation is also performed to check the generality and accuracy of 
this method. 


Furthermore, using the product moments formulae of two correlated standard normal variables 
[Pearson and Young(1918)[, the polynomial equation of p z is derived. For a given p x between two 
correlated random variables, an appropriate value of p z is conveniently determined, and random 
vector with assigned correlation matrix can be generated from standard normal deviates. Through 
numerical examples, it makes evident that non-normal distributions can be well simulated by the 
polynomial model, and random vector with prescribed correlation structure can be generated by 
the proposed procedures. 
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2 Polynomial normal transformation model 

Let X be a continuous random variable with cumulative distribution function (CDF) F(X), let Z 
be a standard normal random variable with CDF 'P(Z). Equate the two CDFs: 



( 1 ) 



X = F~ l [$(Z)\ 

where F _1 (-) denotes the inverse CDF of X. 

Approximate the function F _1 [^(Z)] by a polynomial of degree n: 


( 2 ) 


n 



( 3 ) 


where a,k (k = 0,..., n) are the undetermined coefficients. Setting n = 3 or 5 would yield the TPNT 
or FPNT model. 

Since both F^ 1 (-) and <&(Z ) are monotonic increasing functions, so is the composite function 
F _1 [$(Z)]. According to Weierstrass approximation theorem [Saxe (2002)], the monotonic smooth 
function F -1 [#(Z)] can be well approximated by a polynomial of Z over a closed interval. While, 
the major problem is to determine a set of coefficients (k = 0,... ,n), such that the random 
variable X would be well simulated by the standard normal variable Z. Most often, the researchers 
resort to the raw moment matching method[Fleishman(1978),Headick(2002)], which requires te¬ 
dious mathematical computation. For a polynomial of degree 5, the resulted system of nonlinear 
equations is very complicated (see Appendix A of reference[Headick(2002)]), let alone for a poly¬ 
nomial model of higher degree. In this paper, the PWM matching and percentile matching are 
introduced to simplify the coefficients estimation problem, and the polynomial model is extended 
to higher degree, allowing for the simulation of more probability distributions. 

2.1 Probability weighted moment matching method 
The PWM of X is defined as [Greenwood et al (1979) ]: 


M r ^ s = F{[F(X)] r • X* • [1 - -TOD 

For computational convenience, a particular type of PWM, /3 r = AF. 1 , 0 , is considered: 


( 4 ) 



( 5 ) 


where /(X) is the probability density function (PDF). 

If the PDF of X is unknown, only a set of data is available. Sort the sample into ascending order: 
x\ <•••<£*<•< x m , the unbiased estimate of /3 r is [Hosking et al(1985)Hosking, Wallis, and Wood]: 


m (m — 1 )(to — 2) • • • (m — r) * 


(6) 
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Substitute Eq.(3) into Eq.(5), the PWM of polynomial model is: 

p+oo ( n \ 

Pr = <P r (Z)- [J2 a kZ k \-tp{Z)dZ 

\k =0 / 

n p-\-oo 

= 5 > fc / $ r (Z)-Z k -v(Z)dZ 

fc=o 
n 

= ^2/ a kMr,k, 0 
k—0 

/ +oo 

$ r (Z) ■ Z k ■ ip{Z)dZ 

-oo 


-M,fc, 0 — 


(7) 


where ip(Z) is the PDF. Since d>(Z) and ip(Z) are known functions, M* k0 can be integrated 
numerically. 

For an nth order polynomial model, perform the first (n + 1) PWMs matching, the following 
system of linear equations is established: 


/Mi,0,0 • ■ • M,M ■ ■ • M,n,o\ 


( a °\ 


/AA 


Ml 0,0 • • ■ M; k 0 ■ ■ • M r % i0 


at 

= 

Pr 

(8) 

^Mi,o,o ''' M~ k 0 ■ • ■ M~ n 0 j 


) 


\PnJ 



Calculate the first (n + 1) PWMs of X: 

/ +oo pi 

F r (X) ■ X ■ f(X)dX = / F~ 1 (p)-p r dp (9) 

-oo J 0 

Substitute these PWMs into Eq.(8), solving the system of linear equations gives the coefficients a*,. 

Intuitively, a polynomial of higher degree allows for the control of higher order moments, and 
should give a better approximation of the target distribution. However, as the degree of the poly¬ 
nomial n increases, the coefficient matrix becomes near singular. Below is the determinant of the 
coefficient matrix: 


Table 1: The determinant of coefficient matrix 


n 

Determinant 

n 

Determinant 

i 

0.2821 

7 

3.3321 x 10~ 14 

2 

0.0259 

8 

7.1266 x 10- 18 

3 

8.2291 x 10 -4 

9 

5.7686 x 10 -22 

4 

9.3220 x 10 -6 

10 

1.7762 x 10 -26 

5 

3.8458 x 10 -8 

11 

2.0880 x 10 -31 

6 

5.8600 x 10" 11 

12 

9.4023 x 10" 37 


For a system of linear equations with a near singular coefficient matrix, the solution is very sensitive 
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to small errors in the estimation of k 0 or f} r . Thus, when modeling theoretical distributions, a 
polynomial model of degree higher than 12 is not recommended, otherwise, the near singular coef¬ 
ficient matrix would yield biased values of Ofc, giving rise to an undesirable simulation. Similarly, 
when used for data fitting, a polynomial of lower degree is more preferable, in case of the error 
caused by the high-order PWMs drawn from the samples (In [Zou and Xiao(2014)], a 9th-order 
polynomial model is employed). 


2.2 Percentile matching method 

For a polynomial model, the following requirement should be satisfied: 

x p = a 0 + aiz p H-b a n Zp (10) 

where p is a percentage (p G [0,1]), x p and z p are the corresponding percentiles from the target 
distribution F(X ) and the standard normal distribution <P(Z ): 

x P = F~ 1 (p) z p = <£ _1 (p) (11) 

Select m values of pt , m> (n+ 1), calculate the associated percentiles: z Pi and x Pi (i = 1, • • • , to). 
Based on to pair values of ( z Pi ,x Pi ), the polynomial model can be established by a least square 
method. 

Here, a heuristic method is put forward to establish the polynomial by percentile matching: 

1. Choose 14, 16 and 15 percentage points of pi evenly over the interval [a, 0.01), [0.01,0.99) and 
[0.99,1 — a] respectively, (a is a small value of percentage). 

2. For each percentage p^, calculate the corresponding percentiles of F(x) and <P(z) as stated in 
Eq.(ll). 

3. Using the resulted 45 pair values of ( z Pi ,x Pi ), a 19th order polynomial can be established by a 
least square method. 

To check the performance of the proposed method, the absolute relative error between the 
percentile from the original distribution and the one from polynomial model is calculated: 

x 100[%], x Pi = F~ 1 (pi), 

( 12 ) 

X Pi = GkZ Pi > Z Pi = ^ _1 (Pi)- 
k =0 

Testing for various probability distributions in Matlab, such as: Gamma distribution Gamma(a, b), 
Beta distribution Beta(a,b), F-distribution F(di,d, 2 ), Lognormal distribution lnN(a,b ), Weibull 
distribution Weibull(a,b), Uniform distribution U(a,b), Gumbel distribution Gumble(a,b ), Lo¬ 
gistic distribution Logistic(a,b), T-distribution T(u), Chi squared distribution chi2(k), Rayleigh 
distribution Rayleigh(a) and Exponential distribution Exp(X). the results are shown in Table 2. 

As shown in Table 2, the tested non-normal distributions are well approximated by the 19th 
order polynomial model. Taking the Gamma distribution Gamma(a , b) for example, the error e Pi 
for the family of Gamma distribution with parameters 1 < a < 100, 1 < b < 100 is all less than 
0.92%, demonstrating a satisfactory accuracy. 

However, if the polynomial model is employed for data fitting , a large number of samples should 
be provided to estimate a precise value of the percentile x Pk . With a small to moderate sample size, 
the PWM matching method is more favourable. 


'-Pi 
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Table 2: The absolute relative error of 19th order polynomial model for simulating non-normal 
distributions 


Probability distribution 

Parameters 

Probit range 

[a, 1 — a] 

Maximum value of e Pk (%) 

Gamma 

1 < a < 100 

1 < b < 100 

[ 10- 4 ,1 — 

10 - 4 ] 

0.92 

Beta 

2 < a < 20 

2 < b < 20 

[ 10- 4 ,1 — 

10 - 4 ] 

2.3 x 10~ 4 

Beta 

1 < a < 2 

1 < b < 2 

[ 10- 3 ,1 - 

10" 3 ] 

6 x 10“ 6 

Weibull 

1 < a < 100 

1 < b < 100 

1 

1 

0 

1 

10 - 4 ] 

0.92 

Lognormal 

0.01 < a < 200 
0.01 < b < 4 

[ 10- 4 ,1 - 

10 - 4 ] 

0.060 

F 

4 < di < 100 

4 < d 2 < 100 

[ 10- 4 ,1 — 

10 - 4 ] 

0.090 

Uniform 

a = 0 b = 1 

[10“ 3 ,1 - 

10 - 3 ] 

0.0035 

Gumbel 

<2 = 06=1 

[ 10- 4 ,1 — 

10 - 4 ] 

0.013 

Logistic 

a = 0 6 = 1 

[10“ 4 ,1 - 

10 - 4 ] 

6.6 x 10 -5 

T 

1 < V < 100 

[ 10- 4 ,1 — 

10 - 4 ] 

0.064 

Chi squared 

2 < k < 100 

[ 10- 4 ,1 — 

10 - 4 ] 

0.94 

Rayleigh 

(7=1 

[ 10- 4 ,1 — 

10 - 4 ] 

0.0063 

Exponential 

A = 1 

[10— 4 ,1 - 

10 - 4 ] 

0.95 


3 Evaluating the correlation coefficient 

Suppose X\ and Xi are two correlated random variables. Let X \, X 2 both be approximated by a 
polynomial of degree n: 


Xi — ai : o + aip-Zi + • • • + ai^nZ™ 

X 2 = 02,0 + 02 , 1^2 + • • • + 02 , nZ% 

To ensure a desired correlation coefficient p x between X\ and X 2 , a suitable correlation coefficient p z 

between Z\ and Z 2 should be determined. This problem is known as the Nataf transformation [Liu and Der Kiureghiar 

Lebrun and Dutfoy(2009)] or Normal-To-Anything (NORTA) algorithm[Cario and Nelson(1997), 

Chen(2001)], which requires solving the following integral equation: 

Px (J i (J 2 + P 1 P 2 = E[XiX 2 \ = J J F 1 1 [<P(Zi)]F 2 1 [^{Z 2 )] 4 >{Zi 1 Z 2 , p z )dZidZ 2 , (14) 

where Ui are the mean and standard deviation of Xi (i = 1,2), <f>(Z 1 , Z 2 , p z ) is the joint PDF of 
two normal random variables with correlation coefficient p z . Generally, this equation is difficult to 
be analytically solved. While, using the polynomial normal transformation technique, the problem 
can be greatly simplified. 
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The formulae of the product moments of two correlated standard normal variables are as follows 
[Pearson and Young(1918)]: 


to M/o+M min( ~ s ’ t ) 

E{z} , zr) = &m, y 


(2 Pz) 2 > 


j =0 


(s - j)\(t - j)\{2j)\ 

min(s,t) 


/ „ ^ V . / _ I I t/V I o\ o 

£(z?-+ i z;‘«)= fe (2s+1 2 ) ;|f +1 h £ 


(2pz) 2j 


1=0 


(s-j)!(2j + l)! 


£ (^ S Z^ +1 ) = £ (Z 2s+1 Z|‘) = 0 
(2s)! 


E (Z? s ) = 


i = 1,2. 


2 s • s! 

where s and t are nonnegative integers. 

Then, the following derivation is performed: 

PxO'lO'2 + pL\pL2 = -£(^ 1 X 2 ) = E 




U=0 


iJ=0 


^Ql.o Ol,l 


/ 1 Z 2 • 

Z x Z X Z 2 ••• ZiZJ 


Z r i \ fa 2 ,o\ 


a 2,1 

V a 2,«/ 


(15) 


(16) 


Zj Z‘2 ■ ■ ■ z?z? 

Substitute Eq.(15) into Eq.(16), and perform matrix multiplication, a polynomial equation in p z of 
degree n can be obtained: 

Px&i&2 + P 1 P 2 = b n p " + • + bip\ + • • • + b\p z + bo , (17) 


where 6; (i = 0,..., n) are sums of product of ay*,, and a 2 ,k 2 - A 19th-order polynomial is presented 
in Appendix. 

Then, the correlation coefficient p x is expressed as a polynomial function of p z . As indicated by 
Lemma 3 in reference[Liu and Der Kiureghian(1986)], the feasible correlation coefficient p x cannot 
take any value between —1 and 1, and is located in a subinterval of [—1,1], i.e. —1 < p x < p x < 
1. Based on the derived polynomial equation, the lower and upper bounds of p x can be determined 
by setting p z = —1,1. 

Furthermore, some properties of the function relationship between p z and p x have been proved: 

1. p x is a strictly increasing function of p z (Lemma 1 in reference[Liu and Der Kiureghian(1986)]). 

2. Pz > 0 (< 0) implies p x > 0 (< 0) (Proposition 1 in reference [Cario and Nelson(1997)]), that’s 
to say, p z p x > 0 . 

Thus, for a given correlation coefficient p x between Xi and X 2l the associated p z can be estimated 
by solving the polynomial equation, and a valid solution is restricted by the following conditions: 


— 1 < Pz < 1 and p z p x > 0 


(18) 


As long as Xi and X 2 can be well approximated by the polynomial model, p z for a given p x can 
be estimated with a satisfactory accuracy. 

Suppose X = (Xi, ..., Xi ,..., X m ) T is an m-dimensional random vector with correlation matrix 
Rx, the steps for generating X are as follows: 
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(1) Approximate each random variable Xi by the polynomial model, determine the coefficients by 
the PWM matching method or percentile matching method. 

(2) Calculate p z {i,j ) (* ^ j ) for each entry p x {hj) hr Rx, obtain Rz, which is the equivalent 
correlation matrix in the standard normal space. 

(3) Generate m-dimensional independent standard normal vector: U = (C/i, ■ • • ,Ui , ••• ,U m ) T , 
transform U into the correlated standard normal vector Z = ( Z\ • ■ ■ Zi- • ■ Z m ) T with corre¬ 
lation matrix Rz by performing Z — LU. L is the lower triangular matrix resulted from 
Cholesky decomposition, that is, Rz — LL T . 

(4) Convert Zi into Xi through the marginal transformation Xi = F^ 1 [(I>( y Zi)\ 1 the random vector 
X with the prescribed marginal distributions and correlation matrix Rx is obtained. 

The procedures can be expressed as follows: 


/ Cl \ 

r z =ll t 

( Zl \ 


(Xi\ 


i 




Ui 

Z = LU 

Zi 

X i =Fr 1 (z i ) 

x. t 



\u m ) 


\Zm) 


\x m J 


(19) 


As long as the marginal distribution Fi(Xi) can be well approximate by polynomial model, 
and the correlation matrix Rx is a positive definite symmetric matrix (see Proposition 3 in 
reference[Cario and Nelson(1997)]), the correlated random vector X can be generated by the pro¬ 
posed method. 


4 Numerical example 

4.1 Simulating non-normal distributions 

The polynomial models of degree 3, 5, 11 are employed to simulate the Beta distribution Beta{ 2,2) 
respectively. The coefficients are determined by PWM matching method in Section 2.1, and the 
coefficients of the 11th polynomial model are shown in Table 3. The PDFs of the random numbers 
generated by polynomial models are depicted in Fig. 1. 


Table 3: The coefficients of 11th order polynomial model- Beta(2, 2) 


k 

a k 

k 

a k 

11 

-2.76217093165881 x 10“ 8 

5 

0.00120213182254751 

10 

-8.53583161405554 x lO" 10 

4 

1.41667915745929 x 10“ 7 

9 

1.75657290203781 x lO” 6 

3 

-0.0192416046002625 

8 

1.63734319656302 x 10“ 8 

2 

-6.4707982562276 x 10 -8 

7 

-0.0000555670500494067 

1 

0.265961312977451 

6 

-8.51019206906232 x 10 -8 

0 

0.500000003658777 


10 4 values of pi are evenly chosen in the probit range [0.001,0.999]. Calculate e Pi for each 
percentage pi (i = 1,..., 10 4 ) as stated in Eq.(12), the average value,the minimum value, and the 
maximum value of e Pi are presented in Table 4. 
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Fig. 1: The PDFs of Beta distribution Beta( 2, 2) 


Table 4: The values of e Pi for simulating Beta(2 1 2) 


Degree (n) 

Average(%) 

Minimum(%) 

Maximum(%) 

3 

1.2 

1.8 x 10“ s 

389 

5 

0.15 

8.0 x 10 -6 

54 

11 

4.1 x 10~ 4 

2.0 x 10~ 9 

0.095 


As shown in Fig. 1 and Table 4, due to the ability of controlling higher order moments, a 
polynomial model of higher degree provides a more accurate approximation, the PDF given by the 
11th order polynomial model is in close proximity with Beta distribution. 

The TPNT model cannot provide a good approximation for the Lognormal distribution lnN{ 0,1). 
Here, the polynomial models of degree 5 and 11 are employed to simulate lnN( 0,1). 17 percentage 
points pi are evenly chosen in the interval [0.001,0.999]. Estimate z Pi and x Pi as stated in Eq.(ll). 
With the resulted 17 pair values of ( z Pi ,x Pi ), the coefficients of the polynomial model are deter¬ 
mined by the least square method. The coefficients of the 11th polynomial model are shown in 
Table 5. 


Table 5: The coefficients of 11th order polynomial model- lnN( 0, 1) 


k 

a k 

k 

a k 

ii 

2.70587705587477 x 10~ H 

5 

0.00833335056026873 

10 

3.07095310533251 x 10“ 7 

4 

0.0416665771671574 

9 

2.75245190929274 x 10" 6 

3 

0.166666657138832 

8 

2.46853662677270 x 10" 5 

2 

0.500000016464362 

7 

0.000198405765895619 

1 

1.00000000118744 

6 

0.00138905069923553 

0 

0.999999999545197 


Fig. 2 shows the PDFs of the Lognormal distribution and its polynomial model analogs. Calculate 
e Pi as in the former case, the results are summarized in Table 6. As can be seen, the 11th polynomial 
model gives better performance than the 5th one. 
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Fig. 2: The PDFs of Lognormal distribution lnN(0, 1) 


Table 6: The vaules of e p for simulating lnN( 0,1) 


Degree (n) 

Average (%) 

Minimum(%) 

Maximum (%) 

5 

3.6 

7.2 X 10“ b 

167 

11 

4.6 x 10- 4 

0 

0.087 


4.2 Generating correlated random vector 


Suppose X\ and Xi both follow the Lognormal distribution lnN( 0,1). In this case, the function 
relationship of p x and p z is as follows [Mostafa and Mahmoud(1964)]: 


Px — 


exp{p z ) - 1 
e — 1 


( 20 ) 


For a given p x , p z resulted from Eq.(20) is regarded as benchmark. 

The foregoing 11th order polynomial model is used to simulate X\ and X 2 , of which the coef¬ 
ficients are presented in Table 5. In Appendix, the coefficients a*, (k = 12,..., 19) are set to be 0. 
The results are summarized in Table 7. The function curves are shown in Fig. 3. 


Table 7: The correlation coefficients p z between lnN{ 0,1) and lnN( 0,1) 


px 

p z (Exact value) 

p z (Equation) 

-0.3 

- 0.725 

- 0.725 

-0.1 

- 0.189 

- 0.189 

0.1 

0.159 

0.159 

0.3 

0.416 

0.416 

0.5 

0.620 

0.620 

0.7 

0.790 

0.790 

0.9 

0.935 

0.935 


From Table 7 and Fig. 3, it can be seen that the proposed method yields results with satisfactory 
accuracy. 
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p z 

Fig. 3: The function relationship between p z and p x 


Finally, an example of generating correlated random vector is worked. Suppose X = (Xi, X 2 ,X$) T 
is a three dimensional random vector, X\, X 2 , X 3 follow the standard normal distribution iV(0,1), 
Beta distribution Beta( 2,2) and Lognormal distribution lnN( 0,1), respectively. The desired corre¬ 
lation matrix Rx is: 

/ 1 0.9 0.5\ 

R x = 0.9 1 0.3 

\0.5 0.3 1 ) 

For the normal distribution with mean p and standard deviation a, the coefhcients of the 
polynomial can be regarded as: ao = p, a\ = a and at = 0 (i = 2,..., n). The Beta distribution and 
Lognormal distribution are both simulated by the 11th order polynomial model, the coefhcients are 
tabulated in Table 3 and Table 5 respectively. The equivalent correlation matrix Rz in the normal 
space is: 

/ 1 0.907 0.655\ 

Rz = 0.907 1 0.400 

\0.655 0.400 1 / 

10 6 three-dimensional random vectors are generated as stated in Eq.(19). The correlation matrix 
of samples is: 

f 1 0.900 0.499\ 

R* x = 0.900 1 0.300 

\0.499 0.300 1 j 

Comparing R* x with Rx, it makes evident that the correlation matrix of samples is in close 
proximity to the desired one. Multivariate random variables with prescribed marginal distributions 
and correlation matrix are generated. 


5 Conclusions 

A polynomial normal transformation model is developed in this paper. The PWM matching method 
and the percentile matching method are introduced to determine the coefficients, whereby the 
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polynomial normal transformation model is extended to a higher degree. Comparing to the extant 
TPNT and FPNT models, the proposed model can preserve more statistical information of the 
target distribution and gives more accurate approximation. 

For two correlated random variables, the polynomial equations for estimating p z are also derived. 
It is shown that the original correlation coefficient p x can be expressed as a polynomial function of 
p z . For a given p x , the associated p z is estimated by solving the polynomial equation. The numerical 
examples demonstrate that non-normal distributions can be well approximated by the polynomial 
model, and random vector with prescribed correlation matrix can be generated by the proposed 
method. 
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Appendix 

&19 = 121645100408832000(21 5 19^2,19 
b 18 =640237370572800001,1802,18 
&17 =355687428096000(ai,i7 -t-1 1 loi,ig)(o 2 ,i 7 4* 171o 2 ,i 9 ) 
b\Q =20922789888000(ai,i6 4“ 153oi,i 8 )(o 2 ,i6 4“ 153o 2 ,i 8 ) 

615 =1307674368000(oi,i5 4* 136oi,i7 4- 11628oi,i9)(o2,i5 4- 13602,17 4* 11628o2,ig) 

&14 =87178291200(ai,i4 4- 120o 4 ,i6 4* 9180 oi,i 8 )(o 2 ,i 4 4* 120o 2 ,i6 4- 9180o2,is) 

&13 =6227020800(ai,i3 4- 105oi,i5 4- 7140oi,i7 4- 406980oi,ig) 

(a 2 ,i3 4- 10502,15 4" 7140o2,i7 4- 406980o2,ig) 
b \2 =479001600(ai,i2 4- 91oi,i4 4- 5460oi,ig 4- 278460oi,i8) 

(®2,i2 4- 91a2,i4 4- 5460o2,i6 4- 27846002,is) 
bn =39916800(oi,n 4- 78oi,i3 4- 4095oi,i5 4- 185640oi,i7 4* 7936110oi,ig) 

(a 2 ,n 4* 7802,13 4* 4095o2,i5 4" 185640a2,i7 4" 7936110a2,ig) 

&io =3628800(oi,io 4- 6601,12 4- 3003oi,i4 4- 120120oi,i6 — 4594590oi,i8) 

(a 2 ,io + 6602,12 4" 3003a2,i4 4- 120120a2,i6 4* 4594590o2,is) 
b g =362880(oi, 9 4- 55oi, n 4- 2145 oi,i 3 4- 75075oi,i 5 4- 2552550oi,i 7 + 87297210oi, 19 ) 

(o 2 ,g 4” 5502 ,11 4” 2145a2,i3 4- 75075o2,i5 4- 2552550o2,i7 4- 87297210o2,ig) 
b 8 =40320(oi, 8 + 45oi,io 4- 1485oi, 12 + 45045oi,i 4 + 1351350oi,i 6 4- 41351310oi,i 8 ) 

(o 2 , 8 4“ 45 o 2 ,io 4~ 1485a2,i2 4- 45045a2,i4 4- 1351350a2,i6 4- 41351310a2,i8) 
b 7 =5040(oi,7 + 36oi,g + 990oi,h + 25740a M3 + 675675oi,i 5 4- 18378360ai,i 7 + 523783260ai,ig) 

( 02,7 4~ 3602,9 4“ 990o 2 ,h 4~ 25740o2,i3 4- 675675o2,i5 4* 18378360o 2 ,i7 4- 523783260o2,i9) 
b 6 =720(oi, 6 4- 28oi, 8 4- 630ai,i O + 13860ai,i 2 + 315315ai,i 4 + 7567560ai,i 6 4- 192972780ai,i 8 ) 

(o 2 , 6 + 28a 2 , 8 4- 630a 2 ,io + 13860a 2 ,i2 + 315315a 2 ,i 4 4- 7567560a 2 ,i 6 + 192972780a 2 ,i 8 ) 
b 5 =120(oi , 5 4- 21oi , 7 4- 378ai, 9 + 6930oi,h 4- 135135ai,i 3 + 2837835oi,i 5 + 64324260a M7 4- 1571349780ai, 19 ) 
(a 2 , 5 + 21a 2 ,7 4- 378a 2 ,9 + 6930 o 2 ,h 4- 135135a 2 ,i 3 + 2837835 o 2 ,i 5 + 64324260o 2 ,i 7 + 1571349780a 2 ,i9) 
b 4 =24(ai, 4 4- 15 oi, 6 + 210ai, 8 + 3150ai,i O + 51975 oi,i 2 + 945945ai, 44 + 18918900oi,i 6 + 413513100oi,i 8 ) 

(a 2>4 + 15a 2 ,6 4- 210a 2 ,8 4- 3150o 2 ,io + 51975o 2 ,i 2 4- 945945a 2 ,u 4- 18918900a 2 ,i 6 + 413513100a 2 ,is) 
b 3 =6(01,3 + 10ai, 5 + 105ai,7 + 1260ai,g + 17325ai,n + 270270oi,i 3 4- 4729725ai,i 5 + 91891800ai,i 7 + 
1964187225ai,i 9 )(a 2 ,3 4- 10a 2 , 5 4- 105a 2 , 7 4- 1260a 2 , 9 4- 17325a 2 ,n 4- 270270a 2 ,i 3 + 4729725a 2 ,i54- 
91891800a 2 ,i7 4- 1964187225o 2 , 19) 

5 2 =2(01,2 4~ 6oi, 4 4~ 45ai,6 4* 420oi, 8 4- 4725oi,io 4- 62370oi,i2 4- 945945oi,i 4 4- 16216200oi,ig4- 
310134825oi,i 8 ) (o 2 ,2 4* 6 o 2 , 4 4- 45o 2 ,6 4- 420o 2 , 8 4* 4725o 2 ,io 4- 62370o 2 ,i 2 4- 945945o 2 ,i 4 4- 
16216200a 2 ,ie 4- 310134825a 2 ,is) 

b 3 =( 01,1 4~ 3 oi ,3 4" 15oi, 5 4“ 105oi,7 4- 945oi,g 4- 10395oi,n 4- 135135oi,i3 4* 2027025oi,i54- 
34459425oi,i7 4- 654729075ai,i 9 )(a 2 ,i + 3a 2 , 3 4- 15a 2 , 5 4- 105a 2 , 7 4- 945a 2 , 9 4- 10395a 2 ,u+ 

135135o 2 ,i3 + 2027025a 2 ,i5 4- 34459425a 2 ,i 7 + 654729075a 2 ,i 9 ) 

5q = (oi.o 4~ oi, 2 4~ 3ai, 4 4~ 15 oi,6 4* 105oi, 8 4- 945oi,io 4- 10395oi,i 2 4- 135135oi,i 4 4- 2027025oi,i64- 
34459425ai,i 8 )(o2,o 4- a 2 ,2 4- 3o 2 , 4 4- 15a 2 ,6 4- 105a 2 , 8 4- 945a 2 ,io 4- 10395o 2 ,i 2 4- 135135a 2 ,i 4 4- 
2027025a 2 ,i6 + 3445942502,is) 


(21) 
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